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Abstract 

The phase behavior of hard superballs is examined using molecular dynamics within a deformable 
periodic simulation box. A superball's interior is defined by the inequality \x\'^'^ + \y\'^'^ + [-zl^"^ < 1, 
which provides a versatile family of convex particles {q > 0.5) with cube-like and octahedron- 
like shapes as well as concave particles {q < 0.5) with octahedron-like shapes. Here, we consider 
the convex case with a deformation parameter q between the sphere point {q = 1) and the cube 
{q = oo). We find that the asphericity plays a significant role in the extent of cubatic ordering 
of both the liquid and crystal phases. Calculation of the first few virial coefficients shows that 
superballs that are visually similar to cubes can have low-density equations of state closer to 
spheres than to cubes. Dense liquids of superballs display cubatic orientational order that extends 
over several particle lengths only for large q. Along the ordered, high-density equation of state, 
superballs with 1 < q < 3 exhibit clear evidence of a phase transition from a crystal state to 
a state with reduced long-ranged orientational order upon the reduction of density. For q > 3, 
long-ranged orientational order persists until the melting transition. The width of coexistence 
region between the liquid and ordered, high-density phase decreases with q up to q = 4.0. The 
structures of the high-density phases are examined using certain order parameters, distribution 
functions, and orientational correlation functions. We also find that a fixed simulation cell induces 
artificial phase transitions that are out of equilibrium. Current fabrication techniques allow for the 
synthesis of colloidal superballs, and thus the phase behavior of such systems can be investigated 
experimentally. 
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I. INTRODUCTION 



As the ability to control size, shape, and structure of nanoparticles [l|-l3| and colloids j^ 
improves, computer simulation and theory of hard-particle systems becomes increasingly 
important to the identification of technologically useful properties. Hard convex particles 
have been used as models for simple atomic liquids and solids as a means to connect the 
particle shape and excluded volume to the entropy and to the equilibrium phase diagram of 
a system. The hard-sphere model has a rich history and continues to provide deep insights 
into fundamental physical phenomena |7H15|. However, nonspherical hard particles exhibit 



more complex phase behavior than hard spheres, since the possibility of anisotropic phases 
arises, including smectic, nematic, columnar, and cubatic liquid crystals |l6|. 

The cubatic phase has garnered recent attention because, unlike other liquid- crystalline 
phases, it is characterized by ordering in three mutually perpendicular directions while 
the particles retain translational mobility 17|. Such unusual ordering may lead to novel 
optical, rheological, or transport properties. The cubatic phase has been discovered in 
several hard-particle systems. Cut hard spheres of certain aspect ratios form small stacks 
that a.,g. perpe„d,a,.arly to neighboring staclcs Q. The Onsager eros. a parfe.e eo^sfng 
of three thin rods aligned along orthogonal axes and intersecting at their rnidpoints [19[ and 
tetrapods, hard bodies formed by four rods connected at tetrahedral angles 20| , are examples 
of nonconvex particles that exhibit cubatic ordering. The cubatic phase also arises in systems 
of perfect tetragonal parallelepipeds [2l| as well systems of cuboids, nonconvex particles 
consisting of an array of tangent hard spheres that approximate tetragonal parallelepipeds 
171 . |22| . Perfect tetragonal parallelepipeds have sharp corners and flat faces while the 
cuboid is "bumpy" to approximate friction. Monte Carlo simulation studies of these particles 



revealed a cubatic phase, or parquet phase 



between the liquid and crystal phases 



Pi 
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or aspects ratio other than 1:1:1, that arises 



22|. 



In this paper, we use molecular dynamics (MD) to investigate the equilibrium phase 
behavior and the onset of cubatic ordering in systems of cube-like superballs. A superball 
is a centrally symmetric particle defined by 23 1 

(1) 



Ixl'" + |y|'^ + 1^1'" < 1, 



where x,y and z are Cartesian coordinates and q is the deformation parameter 2J]. Su- 
perballs can take on concave shapes ranging from a cross [q = 0) to the convex octahedron 



(a) 9 =1-0 (b)g=1.5 (c) g = 2.0 (d) 9 = 4.0 (e) g = cx) 

FIG. 1: (Color online) Superballs for certain deformation parameter q. The lines represent the 
three equivalent principal axes. 

(g = 1/2) to a sphere (g = 1) and finally to a cube (g = 00). We focus on the "cube-Uke" 
range 1 < g < 00 to examine the interpolation from spheres to cubes and reveal the role of 
particle shape and curvature on cubatic ordering. As g is increased from unity, the particle 
takes on more cube-like characteristics as edges and corners sharpen while faces flatten. 
Figure [T] displays several superballs with their principal axes for the deformation from unity 
to infinity. 

For a system of superballs, the packing fraction is the fraction of space occupied by the 
particles, = pvsb, where p is the number density, Vsb is the volume of a superball. 



Vsb 



1b(1.?^1±1]b(1' + ' 



2g' 2g 



2g g 



(2) 



B{x,y) = r{x)r(y)/r{x + y), and r(x) is the Euler gamma function 23|. In this study, 
we do not consider the family of superballs for which g < 1. Superballs are well-suited to 
study the effect of the curvature of edges and corners. Although experimentalists have the 
ability to control shape and symmetry of colloidal particles, controlling the curvature may 
be a larger challenge. Scanning electron micrographs of nanoparticles reveal that edges and 
corners are not necessarily sharp like perfect hard polyhedra [3]. Therefore, understanding 
the effects of curvature of hard particles on the phase diagram is not only of fundamental 
interest, but also of practical importance. 

In this study, we detail the phase diagram of hard superball systems as a function of 
g and using molecular dynamics. To characterize the equilibrium phase behavior, it is 
useful to start a system in an equilibrium state and allow the particles to grow or contract 
(or equivalently, allow the simulation cell to shrink or expand). The ideal gas is a suitable 
initial condition to study low-density phases. Fast growth rates applied to liquids of hard 



superballs generate nonequilibirum, randomly jammed packings with novel characteristics 
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26| . However, slow growth rates can also result in nonequilibrium glasses and defective 



crystals and can limit access to complex, high-density equilibrium phases. 

The densest-packing arrangement is a suitable initial condition for high-density phases 
since this arrangement minimizes free energy for a hard-particle system. This study was 



motivated by and made possib 



entire family of superballs 23 



e by the recent development of optimal packings for the 



271. Jiao et al. determined the densest known, and likely 



optimal, packings of superballs |23| and superdisks, the two-dimensional analog |28|. The 



optimal packings of spheres was proved rigorously only recently 



and advances in the 



densest-known packings of aspherical particles, including ellipsoids [30| and the Platonic and 
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32| . will allow researchers to explore the entire phase diagram of 



Archimedean solids 
these hard particles. 

The low-density and liquid crystalline phases of many hard- par ticle systems have re- 



ceived significant research attention including those for cylinders 33 1, sphero cylinders |34| . 
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and parallel superel- 



ellipsoids |35l-l38| , cut spheres [l8|, l39| , tetragonal parallelepipeds [/ 
lipsoids, a perturbation from ellipsoids to cylinders 40|. However, for many of these shapes, 
the optimal packings are yet to be identified {e.g. ellipsoids [30|), and therefore, exploring 
the high-density equilibrium phases remains challenging. 

In the following, we examine the liquid equation of state (EOS), structure of the liquid 
phase, and the virial expansion. For ordered, high-density phases of superballs, we examine 
the crystal branch EOS and use suitable order parameters and correlation functions to 
characterize the entropy-driven phase transitions and cubatic ordering. We find that 

• Excluded volume effects are dominated by edges and corners (Sec. IIVAI) . 

• There exists a phase transition along the ordered, high-density branch of the EOS that 
is associated with changes in long-ranged orientational order (Sec. IIVBI) . 

• The extent of orientational order increases with q at all densities (Sec. |V]), and 

• Fixed system boundaries can produce apparent phase transitions in the ordered, high- 
density systems. (Appendix R|) 

The remainder of this paper is organized as follows. We introduce the deformable box 
molecular dynamics methodology, order parameters, and correlation function in Sec. |TT] and 



review the equations of state for the hard-sphere and hard-cube systems in Sec. IIIII In Sec. 
IIVAI we compare the virial expansion and approximate equations of state to simulation 
data. In Sec. lIVBt we present the EOS for the ordered, high-density and crystal phases 
as generated with deforming box MD simulations and show the onset of cubatic ordering 
as a function of q. The freezing transitions are examined in Sec. IIVCI while the phase 
diagram is illustrated in Sec. IIVDI The structural characteristics are investigated in Sec. 
IVl and discussions are provided in Sec. IVTl In Appendix \^ we discuss the role of system 
geometry on the isotropy of the internal stresses, particularly how fixed boundaries can 
induce apparent, misleading phase transitions. 

II. METHODS 

A. Molecular Dynamics 

Simulation studies of the phase behavior of nonspherical particles typically use Monte 
Carlo simulation methods, primarily due to the simplicity of implementation. One draw- 
back of Monte Carlo methods is the difficulty in achieving collective motion among the 
particles, which is often an important characteristic in systems of anisotropic phases. To 
generate the pressure equations of state for cube-like superballs, we use the Donev-Torquato- 



-|43|. The DTS algorithm generalizes the 



44| to nonspherical, convex particles in- 



Stillinger (DTS) molecular dynamics algorithm [41 
Lubachevsky-Stillinger sphere-packing algorithm ^ 
eluding superballs and elhpsoids. In this algorithm, particles are allowed to grow at a speci- 
fied nondimensional rate 7 (or contract for 7 < 0), or equivalently, the system is compressed 
(or expanded). Contact between particles is predicted using generalized overlap potentials. 
The algorithm allows for shape deformation of the boundary using an approach similar to 
Parrinello-Rahman MD J45l]. In Parrinello- Rahman MD, the "coordinates" of the simulation 
cell are continuously driven by the internal stresses of system which are directly related to 
particle interactions. However, in an event-driven simulation of hard particles, particles only 
interact upon contact and cannot directly interact with the cell. In this paper, a Parrinello- 
Rahman-like algorithm is employed where the "velocities" of the lattice vectors are updated 



after a certain number of collisions based upon the anisotropy of stress tensor 4l|, |43| . The 
mass assigned to the cell is equal to that of the total mass of the particles. 



In the work presented here, we have verified that the pressure tensor remains isotropic 
on average when employing the deforming box algorithm. Although this algorithm may not 
rigorously sample an isostress or constant-pressure ensemble, it is a reasonable approxima- 
tion. We note that simulations using fixed boundaries exhibited pressure tensors that were 
anisotropic and gave rise to nonequilibrium phase behavior as detailed in the Appendix. 

We limit our study to g < 4 since the algorithm is numerically unstable for q > 4.0 



431]. Periodic boundary conditions were employed. The reduced pressure is defined as 
Z = p/ pksT, where p is the system's pressure, p is the number density, and ksT is the 
usual energy scale for hard-particle systems. Particles are of unit "diameter," the surface- 
to-surface distance of a chord along one principal axis {i.e. the shortest chord). 

To obtain the liquid EOS, particles were placed randomly in a low-density configuration 
inside a cubic box. They were given random linear and angular velocities and allowed to 
grow at a specified rate 7 until the system reached a defined pressure. For the crystal branch, 
particles were initialized in a slightly expanded form of the densest lattice configuration, with 
the number of particles A^ chosen to be commensurate with the lattice, assigned random 
linear and angular velocities, and simulated using a contraction rate 7 < 0. 

The densest-known packings of cube-like superballs occur in one of two families of Bravais 
lattices, denoted as Co and Ci [23]. For 1 < g < 1.1509, the densest packings of superballs 
are achieved with the Cq lattice, a perturbation of the FCC lattice. For superballs with 
q > 1.1509, the Ci lattice, a deformation of the simple cubic lattice, represents the densest 
arrangement. Since the MD algorithm is slow at high densities due to the high frequency of 
particle coUisions, the initial crystal configurations were unsaturated, typically near 80% of 
the maximum possible packing fraction. 

Growth rates in the range 10~^ < I7I < 10~^ were utilized. The simulation data for 
spheres was compared to widely-accepted data and we find that equilibrium was well ap- 
proximated with I7I < 10^^. Obtaining a full sweep of the density at I7I = 10~^ required 
over two weeks of computation time with 1000 particles, and therefore significantly slower 
growth rates over the entire density range were not practical. 

We more closely examined parts of the phase diagram in which phase transitions were 
evident by running the algorithm with slower rates using near-equilibrium configurations 
as the initial conditions. In some cases, we averaged over const ant- density (7 = 0) MD 
trajectories consisting of nearly 10^ collisions per particle. We find that the phase transitions 



in the 7 = cases are slightly sharper than those in cases using slow growth or contraction 
rates but generally occur at the same densities. Growth rates of 7 > 10~^ us ually produced 
jammed, metastable structures. These are explored in separate studies 
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261. For each 



g, we obtained multiple, independent sweeps of the density for both the liquid and solid 
branches to examine the variability of the results. The results presented in this paper were 
obtained using 1000 particles. We have varied the system size between A^ = 216 to 1728. 

It is important to point out that no simulation method can effectively determine rare 
events. Even detailed free-energy calculations may encounter the inability to find states as- 
sociated with rare events. Unfortunately, hard-particle MD requires serial calculations and 
advances in the parallelization of these algorithms are required to produce longer trajecto- 
ries. Until then, these algorithms are the most efficient use of computer resources and can 
complement Monte Carlo methods that lack the ability to capture dynamics of collective 
motion. 

B. Quantifying Order 

Superballs have three equivalent principal, mutually orthogonal axes labeled A, B and 
C. For each axis j = A, B, or C of superball i, there is an associated unit vector 
Ujj = [uij^x,Uij^y,Uij^z]- For a nematic-forming system, there is at least one "director," 
n = [nx,ny,nz], which represents the most aligned direction in a system. For sufficiently 
anisotropic (large q) superballs, one might expect systems to have at least one nematic 
direction in a low-density phase and three orthogonal directions in a crystal phase. 

Order parameters are useful to characterize the local and global order in a system of 



particles and several have been employed for particles with cubic symmetry 2l|. We find 
that the nematic and cubatic order parameters, S2J and 5*4 respectively, are the most useful 
scalar metrics for quantifying orientational order. For cubatic ordering, a nematic order 
parameter can describe ordering in each direction A, B and C. The nematic order parameter 
for a particular set of axes j, S2J is defined as 



^Ef|K-^^l'-^) (3) 



5*2 ,• = max , ^ . , , „,,, --, , 

i ^ 

where A^ is the number of particles, Ujj is a set of particle axes, and n^ is the director for 
direction j. The solution to Equation ([3]) can be found by solving the eigenvalue problem 

8 



Arij = Arij where 

i 

where 6 is the Kronecker deha function. S2J is the maximum eigenvalue Xmax and the 
nematic director vector rij is the eigenvector associated with Xmax- 

Since all principal axes of a superball are equivalent, we must "relabel" the particle 
axes prior to calculating 5*2^- in order to obtain meaningful data. A set of three mutually 
orthogonal unit vectors is chosen as a reference, typically the principal axes of one randomly 
chosen particle or the standard laboratory axes. For each particle, we relabel the particle 
axes based on the best alignment with the reference system. For example, we identify the 
axis of each superball that is best aligned with the [1, 0, 0] vector, label these axes as A axes, 
then continue with the [0,1,0] vector and B axes. The remaining axes are labeled as C 



axes. A schematic of this procedure is shown in Ref. [ITJ]. The relabeling scheme introduces 
artificial correlations so that S2J in an isotropic system is approximately 0.55. For perfect 
cubatic ordering, S2J = 1 in each of the three orthogonal directions. Here, we report 5*2 as 
the maximum of the S2,/s since the >S'2,j's were nearly equal to each other in all of the cases 
considered. This suggests that we encountered isotropic and cubatic phases and not phases 
with strict uniaxial or biaxial order. 

The cubatic order parameter 6*4 is a more appropriate scalar metric for ordering in three 
orthogonal directions and is defined as 

5*4 = max — j-TT y^ (35|ujj ■ n|'' — 30|ujj ■ n|^ + 3) . (5) 

Here, n is a unit director for which 5*4 is maximized. The prefactor of 1/14A arises from 
the accounting for the 3 A principal axes and normalizing 5*4 to unity for perfect alignment. 
This can be formulated into a eigentensor problem as in Ref. [39]]. However, we use an 
approximate solution. We choose the maximum 5*4 from a large set of trial directors n. 
Here, we take the set of trial directors to be the set of all particle axes Ujj, providing 3 A 
trial directors. We report the maximum 5*4 from this set of trial directors. For perfect 
cubatic ordering, 5*4 = 1, and for a system with no long-range cubatic order, 6*4 = 0. 

We use an orientational correlation function 6*4 (r) to measure the mutual alignment of 
particles as a function of distance r between particle centers. As a specific instance of a 



general class of orientational correlation functions [l8|, 6*4 (r) is defined as 



G,{r) = A (35 [u„,(0) ■ u,,(r)]^ - 30 ^(0) ■ u,,(r)]' + 3> , (6) 

where the < ■ ■ • > denotes the average over all axes j and particle pairs a and b. The 
prefactor of 3/14 accounts for all nine combinations of axes between particle pairs which is 
similar to the normalization of 5*4. In the limit that r — t- 00, G^lr) approaches Si- This 
function allows us to determine local correlations whereas 5*4 determines global order. 

The local alignment between "neighbor" particles can be characterized by the angular 
distribution function f{6) defined for 0° < 6 < 180°. The probability of finding two neighbor 
particles whose axes are aligned at an angle between 6 and 6 + d6 is given by f{9)dd. We 
consider two particles to be neighbors if they are separated by less than 1.35 diameters, 
since for most particles, this is between the first and second neighbor shells of the associated 
crystal. However, f{9) is relatively unaffected by small variations of the cutoff radius. For 
the perfect crystal of cube-like particles, f{0) = ^S{9) + ^6{9 — 90°), where 6{9) is the Dirac 
delta function. At equilibrium, f{6) is symmetric about 6 = 90°. 



The order parameters derived from the spherical harmonics [46[ were calculated but do 
not add significantly to our analysis. The use of a different reference crystal for each q 
precludes our ability to compare superballs with different q. Also, changes in these order 
parameters as a function of (p closely matched changes in 5*2 and 6*4. 

To examine translational order, we use the radial distribution function (72 (r), which is the 
normalized pair density distribution function such that for a disordered system it tends to 
unity for large pair separation r. In addition, one-dimensional particle distribution functions 
were calculated in each of the three associated one-dimensional directions of the crystal. 
When simulating the melting of the crystal, these particle distribution functions remained 
periodic until the crystal melted. This confirms that lower- dimensional translational order, 
such as that found in columnar or smectic phases, was never present. 

C. Virial CoefRcients 

The virial expansion for hard particles in terms of reduced pressure Z is given by 

00 
Z = p/pkT = l + J2 (Bi/v'-') <P\ (7) 



j=2 



10 



where p is the pressure, Bi is the i^^ virial coefficient, and v is the volume of a single 
particle. For hard, convex particles the second virial coefficient is known analytically as B2 = 
RS + V [47| , where R and S are the radius of mean curvature and surface area of a particle, 
respectively. Evaluating analytic expressions for these quantities is highly nontrivial, as 
was the case for ellipsoids [48|], but numerical calculation of the second and higher virial 
coefficients is straightforward given a suitable overlap function. We calculate the ffist few 
virial coefficients using Monte Carlo integration |8|. Trial configurations were generated 



using the method of Ree and Hoover lOj. For B2, 1.5 x 10^ random trial configurations were 
used. For B^, random configurations were generated until 2x10^ configurations satisfied the 
condition that particle 1 overlaps particle 2, 2 overlaps 3, and 3 overlaps 1. For B4, random 
configurations were produced until 2 x 10^ configurations satisfied the condition that particle 
1 overlaps 2, 2 overlaps 3, 3 overlaps 4, and 4 overlaps 1. The standard deviations of ten 
subaverages were less than 0.5% of the virial coefficient. The algorithm was tested against 
known results for hard ellipsoids kol. We have calculated hard-cube virial coefficients to a 

n 

higher accuracy than in Ref. [50| using the separating axis theorem to check for overlaps. 



III. REFERENCE SYSTEMS: HARD SPHERES AND HARD CUBES 

Since superballs interpolate between a sphere and cube for q ranging from unity to infinity, 
we review the phase behavior of these reference systems. The hard-sphere liquid and crystal 
have been well-characterized by simulation and theoretical treatments. The hard-sphere 



EOS as generated by the DTS algorithm with 7 = ±10^ and g = 1 is shown in Fig. 2(a 
Although the DTS algorithm does not achieve true equilibrium, which requires 7 = for long 
MD trajectories, the slow growth rate approximates equilibrium well. The DTS algorithm 
correctly produces the well-known and widely-accepted Carnahan-Starling EOS for liquids 
[12] and the empirically-derived Speedy EOS for the FCC crystal 13|. 



The freezing and melting points of a hard-sphere system are = 0.490 and (f) = 0.545, 



respectively, with a coexistence pressure Z = 11.48 [51i- In Fig. 2(a), one can see that the 
DTS algorithm produces a ffist-order phase transition at = 0.551 to a partially crystalline 
system when starting the system in a random initial configuration and allowing the particles 
to grow slowly at 7 = 10^^. This packing fraction is repeatable over independent runs and 
consistently occurred between (p = 0.545 and 0.553. Halving the growth rate to 7 = 5 x 10~^ 
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FIG. 2: (Color online) (a) Pressure EOS for hard spheres generated by the DTS algorithm with 
N = 1000. The liquid and partially crystalline branches were obtained starting from a random 
initial configuration and 7 = 10~^ while the crystal branch was obtained by initializing the system 
in FCC arrangement and 7 = — 10~^. The coexistence line is that of Refs. 
equation of state and coexistence line for freely rotating hard cubes from Refs. |2ll. l52|. 



.(b) Pressure 



(not shown in figure) resulted in a freezing event at = 0.547. Starting in the FCC crystal 
arrangement and using 7 = —10^^, the simulation data traces the Speedy EOS and the 
system shows a first-order transition at = 0.495. The density at which the phase transitions 
occurred showed little variability, occurring between cj) = 0.496 and 0.498 across multiple 
runs. Halving the rate to 7 = —5 x 10^^ yielded a transition at = 0.496. The algorithm 
evidently is appropriate for determining the densities at which phase transitions occur and 
identifying coexistence regions. The DTS algorithm, however, cannot explicitly identify 
the coexistence pressure, which requires extensive free-energy calculations. An equal-area 
construction of the pressure- volume equation of state may provide an approximation to the 
coexistence pressure. 

There has been considerably less attention devoted the hard-cube system. However, 



several studies have elucidated the EOS, phase transitions, and orientational orderiiig. 



ar- 



allel hard cubes undergo a continuous melting transition from a crystal to a liquid [53|, |54 |. 
The EOS for parallel hard cubes has several gentle curvature changes but no evidence of a 
first-order phase transition |55j. 
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However, we are interested in the EOS and phase transitions in systems of freely-rotating 
hard cubes. John et al. obtained the phase diagram for freely rotating hard cubes by Monte 
Carlo methods 2l| and their EOS is shown in Fig. 2(b) A first-order melting transition 



was first calculated to be between cf) = 0.45 and 0.52 53|. More recent studies reveal a 



narrower coexistence region between an isotropic and a cubatic phase to be in the region 



0.437 < < 0.495 [2l|. More interestingly, John and coworkers suggest that there is 
a cubatic-crystalline transition between (p = 0.634 and 0.674. However, the nature of the 
iransition was not well-characterized and necessitated free-energy calculations for verification 



2l| . The DTS algorithm is numerically unstable for superballs with large q and therefore. 



these Monte Carlo calculations are the best available reference data. 

IV. EQUATIONS OF STATE 

A. Isotropic Liquid Phase 

The low-order virial coefficients describe the low-density EOS of a system. The second, 
third, and fourth virial coefficients, shown in Fig. [HI monotonically increase with the defor- 
mation parameter q. This is expected for B2 and -B3 given that the cluster integrals involved 
in the calculations effectively measure excluded volume, and a larger particle generally yields 
a larger excluded volume. The effects on B4 and higher-order coefficients are not intuitively 
obvious, since these coefficients are not necessarily strictly positive for hard particles. The 
virial coefficients are closely approximated by an exponential equation, 

R 

aiexp{-bi/q) + Ci. (8) 



yi-l 



The values of the parameters obtained by nonlinear least-squares regression are given in 
Table [H These parameters were the best obtained, though they are not guaranteed to be 
optimal due to the nonlinearity of the fit. The insets in Fig. |3]plot 



hi 



- Ci 



yi-l 



(9) 



versus 1/g to show that the scaling of Wi is approximately linear with 1/g, except near the 
sphere point q = 1 where the fit tends to degrade. 

A significant portion of the excluded volume evidently arises from the particles' edges 
and corners. Sharpening of edges and corners increases the excluded volume in an isotropic 
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FIG. 3: The (a) second, (b) third, and (c) fourth virial coefficients for hard superballs fit to the 



equation -j^ = ajexp(— 6j/q) + Cj. The inset shows Wi 
values for the fitting coefficients is provided in the text. 



-fin 



. The table of 



TABLE I: Coefficients of fit for Bi/v'~^ = ai exp{-bi/q) + a. 



B2/V 



B^y 



B^/v^ 



ai 
hi 



1.638 
2.770 
3.869 



8.759 
3.262 
9.564 



24.077 
3.794 
17.715 



phase faster than it increases the actual vokime of a particle. When visually comparing 
a superball with a moderately large q value, say g = 4, to the perfect cube (see Fig. [1]), 
one might surmise that because the two particles have such similar appearances, the two 
particles would have similar behavior in the liquid phase. However, when comparing B2/V, 
the relative excluded volume, for these two particle shapes, the superball with g = 4 is 
closer to a sphere than to a cube. The edges and corners evidently are dominant features 
contributing to excluded volume effects. 

Fig. nil shows the simulation data for several q values for low densities. We compare our 



simulation results to the Nezbeda EOS 
EOS 



56 



571], a modification of the Carnahan-Starling 



12| for convex hard particles. Using only the nonsphericity parameter a = RS/v, the 
Nezbeda EOS is given as 



1 



+ 



Sac 



+ 



3aV - a(6a - 5)(j)^ 



(10) 



1-0 (1-0)^ (1-0)^ 

The Nezbeda curve. Fig. Ht, follows simulation data of superballs along the entire liquid 
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FIG. 4: (Color online) Liquid equation of state for q = 1.5, 2.5, and 3.5 from (a) simulation and (b) 
the Nezbeda EOS. When overlaid, the Nezbeda EOS is accurate for low densities but accelerates 
its divergence at higher densities, typically for (p > 0.25 and for q > 2.5. 

branch for q < 2.5. For q > 2.5, the Nezbeda curve follows the simulation data at low 
densities, slightly underestimates pressures at moderate liquid densities, and slightly over- 
estimates pressures at higher densities. For example, with q = 2.5, the Nezbeda curve is 
accurate for < 0.20, underestimates the pressure for 0.20 < < 0.47, and overestimates 
the pressure for (p > 0.47. Using local polynomial fits to the simulation curves, we can com- 
pare the pressure values of the simulation curves to those of the Nezbeda equation of state. 



For q = 2.5, Zgim — Z^^f.^ is less than 0.059 at = 0.2, while for = 0.35, Z^ 



ZNez is 



less than 0.138. At = 0.48, Zgim — Z^ez is about -0.086. As shown by the simulation data, 
superballs that are visually similar to cubes have pressures midway between that of spheres 
and that of cubes, demonstrating the effects of sharp edges and corners on the low-density 
EOS. 

B. Ordered, High-Density Phases: Melting 

To obtain the ordered, high-density equations of state, systems were initialized in the 
optimal lattice configuration and a contraction rate of 7 = — 10~^ was applied until the 
system entered the fluid phase. Figure [5] illustrates the resulting equations of state for various 
values of q. As seen in the figure for q < 3.0, there exist two apparent phase transitions. 
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FIG. 5: (Color online) Crystal branches and melting transitions for various q values obtained using 
the contraction rate 7 = — 1 x 10^^. The transition from the K phase to Q occurs at a lower 
density for increasing q until it is absorbed entirely by the K phase. 

which result in three distinct phases. We define these phases, going from highest density to 
lowest, as the crystal {K), cubatic (Q), and liquid (L) phases. For q > 3.0, only one phase 
transition is evident from the pressure EOS, separating the K and L phase. 

The K phase is characterized by long-ranged translational and orientational order, quan- 
tifiable by the appropriate order metrics and distribution functions. We characterize the 
Q phase as having a moderate degree of long-ranged orientational order compared to the 
crystal phase. Particles are loosely braced and have an 5*4 value above that of a dense liquid 
but less than that of a crystal, between values of 0.05 and 0.20. In this phase, crystalline 
translational order remains. While some have characterized the cubatic phase as having no 
long-ranged positional order 33|, others also consider those phases with an "intermediate 



degree of translational order" as cubatic phases [17|, |22|. We choose to use the latter defini- 
tion of "cubatic," since this was the definition employed for the use of hard cubes. The L 
phase lacks long-ranged positional and orientational order. 

The nematic and cubatic order parameters associated with the ordered, high-density 
phases are shown in Fig. [61 In all cases, the curve for 5*2 exhibited more noise than the 
corresponding 6*4 curve, possibly due to the "relabelling" scheme. Regardless, the behavior 
of both order metrics follows closely with the pressure equations of state. When the pressure- 
density curve is discontinuous or shows a change in curvature, the order parameters exhibited 
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FIG. 6: (Color online) Nematic and cubatic order parameters 5*2 and 5*4 for q = 1.5, 2.0, 2.5, 3.0 
and 3.5. with 7 = —10^^. The inset shows the behavior at the melting transition for q = 2.0 and 
3.0, though the behavior is similar at this transition for all q. The labels on the inset are the same 
as the larger graph. 

a corresponding change. 

Beginning in the densest crystal and moving down along the EOS, the slope of the reduced 
pressure is discontinuous at the K-Q transition. This transition point initially occurs at high 
densities for small q. For increasing q, this transition occurs at lower densities and by g = 3.5, 
the transition vanishes. As seen by the order parameters in Fig. [6l this transition results 
in a reduction in cubatic ordering. While the S'4-density curve remains continuous at the 
transition, the slope does not, similar to the behavior of the pressure-density curve. Although 
the Q phase has less orientational order than the K phase, there remains significant long- 
ranged cubatic order. At the melting transition from Q to L or i^ to L, the nematic order 
parameters are discontinuous as shown in the inset of Fig. [61 

The K-Q transition was observed in the EOS for 1.3 < q < 3.0, however, we strongly 
suspect that similar transitions occur for smaller q and high pressure. Unfortunately, sim- 
ulating the system for q = 1.2 and (f) > 0.72 proved challenging, both in the stability of 
the code and achieving near-equilibrium behavior. Because of the monotonic behavior of 
the location of the K-Q transition with respect to q, we suspect this to continue for q near 
the sphere point. For q > 3.0, we do not observe a K-Q transition, although it is possible 
that a K-Q phase transition exists but is of a higher order. Higher-order phase transitions 



17 



were suggested for hard cubes [21j, though in the case of nearly cubic superballs, there is no 
discontinuity in the cubatic order parameters aside from the first-order transition associated 
with melting. In addition, we have shown that edges and corners play an important role 
in the liquid EOS, and therefore, it is not unreasonable to expect that the high-density 
behavior of nearly cubic superballs (g = 4.0) can deviate from that of hard cubes. 

As the density is reduced in the Q phase, the cubatic order parameters smoothly decreases. 
At the Q-L transition, the pressure jumps while 5*2 and 5*4 drop close to the values associated 
with random rotations. The Q-L transition is clearly first order and is present for all q 
tested. The density at which melting occurs increases monotonically from cf) = 0.494 for 
hard spheres up to 0.536 for q = 4.0. While reducing the density along the high-density 
equation of state, the translational order appears to drop continuously. The peaks in (72 (r) 
maintain crystal-like characteristics. We discuss translational order in greater detail in the 
next section. Figure [7] shows a representative sample of the system at several densities in the 
K and Q phases. The decreased orientational order associated with the Q phase is not easily 
distinguishable from the visual inspection but is clearly evident by the order parameters. 

The decrease in the density associated with the K-Q phase transition for increasing q is 
attributed to the broadness of the edges and corners of particles as they relate to particle 
rotations. For each q in the K phase, particles are braced against rotations. Particle rota- 
tions require large local fluctuations in the density to allow for particle flips. As the relative 
lattice spacing is increased, the rotational mobility of the particles increases. With increas- 
ing q, particles have sharper corners that are more effective in bracing against rotations for 
larger density ranges. The required local volume necessary for particle rotations increases 
with more cube-like shape, and therefore we observe a larger K phase for increasing q. 

The increased rotational mobility of the Q phase also imparts larger internal stresses 
than in the K phase, as observed by the greater slope of the pressure-density curve in the 
Q phase than the K phase at the K-Q transition. In Sec. IVIl we detail some of the tests we 
performed to ensure that these results were not related to system-size, boundary, or kinetic 
effects. 
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(a) = 0.626 




(b) 4) = 0.600 



(c) (f) = 0.558 

FIG. 7: (Color online) Two views of particle configurations from the crystal branch for q = 2.0 for 
a system containing 512 particles. The dark lines are the boundaries of the simulation box. One 
can observe significant translational order in each. Orientational order is reduced as the density is 
decreased, but translational order is maintained. 

C. Ordered, High-Density Phases: Freezing 

The freezing transitions for superballs were examined by applying a slow growth rate 
to particles initialized in a low-density liquid. Figure [8] shows the pressure as a function 
of packing fraction for two values of q along with the crystal branch equation of state. In 
the figure, two growth rates are displayed, 7 = 10~^ and 5 x 10^^. We find that upon 
freezing most systems order into a partially crystalline structure representative of a L-Q 
transition. Across this transition, the pressure drops while the order parameters increase 
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discontinuously. At high pressures, the particles increase the face-to-face contacts, but the 
nematic directors of the separate crystal regions destructively interfere. Therefore, the order 
parameters cannot achieve values as large as the equilibrium branch at the same density. The 
pressure of the partially crystalline phase is greater than that of the equilibrium structure, 
which is evident in Fig. [HI These systems have characteristics of translational order of a 
crystal but have a large grain boundary or vacancies. At very high pressures, the grains are 
eliminated by the shear deformation of the box. The structural characteristics of jammed 



superballs are detailed in Ref. [26 1. 



The Q-K transition is sensitive to the extent of crystallization occurring at the L-Q 
transition. In several cases, the system slightly overshoots the Q-K transition while in other 
cases, the system never shows signs of a Q-K transition. One might expect that using slower 



growth rates would alleviate this phenomenon. However, as seen in Fig. 8(a) , slower growth 
rates are not necessarily better at achieving the Q-K transition and tracing the crystal 
branch EOS. The density at which the crystallization from the L phase to the Q phase 
occurs is relatively insensitive to the growth rate for the slow rates used in this study. It is 



of note that substantially faster growth rates produce amorphous, jammed structures [26|. 
For q > 2.5, we find that the initial freezing from a liquid phase occurs at a density 
higher than that of the Q-L transition along the ordered, high-density branch which is 



shown in Fig. 8(b) In these cases, both growth rates, 10~ and 5 x 10~ , resulted in a 
second transition, presumably a Q-K transition, and then approximated the curvature of 
the crystal branch. The pressure remained higher than that of the perfect crystal. The order 
parameters increase continuously after this freezing transition. For q > 3.5, we find that 
the L-K transitions along the growth and contraction branches approach one another. It 
is possible that relaxation times are faster for these q values and that transition to a liquid 
state is more easily achieved. 

D. Phase diagram 

As with the case of spheres, we consider the density region between the L-Q transition 
of the growth branch and the Q-L transition of the contraction branch to be coexistence 
between these two phases. In Fig. |H1 this coexistence region is presumed to be the region 
between sharp phase transitions along the growth and contraction branches. Using the 
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FIG. 8: (Color online) Freezing transition for q = 1.5 and 2.5 and 7 = 10^^ (blue), 5 x 10~^ (green), 
and —10^^ (red). In all cases tested, a freezing transition resulted in crystals with vacancies and/or 
grain boundaries. 

results of the DTS algorithm and the quantification of orientational and translational order, 
the approximate phase diagram is shown in Fig. [9l 

The black circles represent the densities of the Co and Ci lattices at maximum packing. 
The blue diamonds represent the density of the K-Q transition as found along the crystal 
branch. The precise boundary was determined using a linear regression of data points on 
either side of the transition and calculating their intersection. The red squares represent 
the first-order phase transition associated with the melting of the crystal. It is the last 
density at which the crystal was stable in our simulations. The green diamonds represent 
the freezing transition by allowing particles to grow slowly. The data points along the cube 
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line (1/g = 0) represents the transitions identified by the authors of Ref. 2l| with the brown 
triangle being the crystal-cubatic transition they identified. 

While the K-Q transition is nionotonic in g, the L-Q transition is not. For q > 2.5, 
the K-Q transition of the contraction branch lies at a density below the L-K transition on 



the growth branch as in Fig. 8(b) , The apparent tie line grows slightly as q increases from 
unity and begins to narrow significantly around q = 2.5. By q = 4.0, the K-L transition 
approaches the same density as the L-K transition. There are several possible explanations 
to the narrowing of the apparent coexistence region near q = 4.0. We may not observe 
the possible coexistence region due to the relaxation times in the system. The free-energy 
barrier associated with melting may be reduced for increasing q and thus the crystal may be 
difficult to stabilize in a coexistence region. Alternatively, we may be observing microphase 
separation where the melting into a liquid phase is preferable to maintain the crystal phase. 
The system may also be large enough where certain domains of the system are liquid while 
others remain crystal, although we do not directly observe evidence for this phenomenon. 
Further refinement of the phase boundary would be necessary to answer these questions. 

Detailed free-energy calculations are necessary to determine the precise phase bound- 
aries, refine the approximation of the phase diagram shown in Fig. [9l and determine the 
coexistence pressures. One could approximate the coexistence pressure using an equal area 
construction of the two pressure branches. We expect that the coexistence pressures would 
be nonmonotonic in q as well. The coexistence pressure for the sphere is Z = 11.48, and for 



the cube, Z ~ 12.5. Using Fig. 8(b) as an example after rescaling the axes appropriately, 
one can estimate the coexistence pressure to be much higher than 12.5. Ultimately, the 
free-energy calculations would definitively determine such a conjecture. 

V. STRUCTURAL CHARACTERISTICS 

A. Liquid Phase 

The liquid phase of superballs lacks long-range orientational and translational order but 
can show significant local order. Comparing systems at the same but for various values 
of q, we find that increasing q yields increased local orientational order but diminished 
local translational order. The radial distribution function for dense liquids at various q for 
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FIG. 9: (Color online) Approximate phase diagram for superballs spanning from cubes to spheres. 

the common (f) = 0.54 is shown in Fig. [TOT a). The first peak in the (72 (r) is reduced with 
increasing q, hkely due to the interaction anisotropy building up in the hquid. Because of 
the anisotropy of the particle, contacting neighbors are not restricted to r = 1 as they are 
for spheres, and thus the first peak in g2{r) is lower than that for spheres. 

The strong peak in G^lr), defined by Eq. [6] at r = 1 for all particle shapes, shown by 
G^lr) in Fig. fTOT b). shows that particles with cube-like shape have a strong preference to 
align orthogonally with contacting particles. Evidently, the existence of cube-like shape 
is sufficient to produce preferential cubatic alignment at contact, even for particles with 
q just above unity and whose cube-like shape is difficult to see visually. The degree of 
curvature determines the range of order in the liquid phase. For larger q, the orientational 
correlations persist for up to three diameters. Sharper corners effectively "brace" particles 
against rotations. The angular distribution function f{6), Fig. [TOJfc). reinforces the notion 
of bracing. With sharper edges, nearest neighbors have a greater preference to align along 
mutually orthogonal directions in the liquid phase. 

B. High-Density Phases 



While 5*4 measures the global orientational order in the system, the orientational corre- 
lation function 6*4 (r) provides insight in the locality of cubatic ordering. Figure [TT] shows 
G^lr) for several q values. For q = 1.5 in the K phase. Fig. 11(a) and = 0.68, the long- 
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FIG. 10: (Color online) (a) Radial distribution function g2{T)-, (b) orientation correlation function 
Gii{r), and (c) angular distribution function f{6) for various q in the liquid phase at (/> = 0.54 for 
N = 1000. The peak of the radial distribution function becomes smaller with increasing q while 
the first neighbor peak moves farther due to excluded volume effects. The local alignment increases 
with increasing g, but diminishes within four particle diameters. Mutual alignment of neighboring 
particles increases sharply for increasing q. 

ranged orientational order is evident by the large value of 6*4 (r) for nearly all r. Particles 
with "face-to-face" contacts at r = 1 and those at pair distances associated with lattice sites 
are more strongly aligned than those that deviate from lattice sites. However, those pairs of 
particles that are separated by 1.35 diameters are misaligned, evidenced by 6*4 (r) which is 
nearly zero at this distance. Once the system enters the Q phase, orientational correlations 
remain local, although 5*4 suggests there remains some global cubatic order. These local 
orientational correlations are slightly longer ranged than in the liquid phase, which results 
in a larger value of 5*4. 

For larger g, we observe similar trends, though particles are better stabilized at all pair 
distances. The bracing of particles with sharper edges prevent rotations at the pair distance 



of 1.35 diameters. For q = 2.5, Fig. 11(b) shows that 6*4 (r) has an initial maximum and 
minimum, followed by weak maxima and minima associated with neighbor distances. The 
general shape of the G4 curve is maintained through all densities associated with the K, aside 
from general loss of long-ranged order. In going from the K to Q phase, the orientational 
correlations at long-ranged nearly vanish in the same manner as that of g = 1.5. Fig. 



11(c) shows that the orientational correlations for q = 3.5 are similar to those at g = 2.5. 
Surprisingly, the contact value of G^^r) was relatively unaffected by the density in all cases. 
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FIG. 11: (Color online) Orientational correlation function G4(r) for various (p. 

The degree of curvature and density have little influence on the ordering at contact, but 
rather it is the presence of cubic symmetry that produces alignment at r = 1. 

To further examine the local ordering environment, we used the angular distribution 
function f{6), shown in Fig. [12] for several q. This measures the correlations between the 
alignment of principal axes for neighboring particles. Local cubatic ordering is present when 
the function has three maxima but has minimal cubatic ordering when there is a single 
maximum. When initialized, the crystals are perfectly aligned and the angular distribution 
function is /(0°) = 1/3 and /(90°) = 2/3. Since the axes are labelled [i.e. all of the A 
axes are in the same direction), the crystal system will have an f{6) that is asymmetric in 
6 at high densities. At equilibrium, particles will have flipped sufficiently to mix the labels 
and f{6) will become symmetric about 90°. The data shown in Fig. [12] illustrates that this 
symmetry is easily achieved at these densities even in the K phase. 



For q = 1.5, Fig. 12(a) shows that superballs in the K phase maintain strong neighbor 



bracing. Orthogonal alignment among neighbors is the most probable alignment. Entering 
the Q phase reduces the peak at 6* = 90° and eliminates the other local maxima. The shape 
of the curve is similar to that of the dense liquid phase. 



For q = 2.5, Fig. 12(b) demonstrates that the K phase prevents nearly all neighboring 
particle axes from aligning at 45° since the depth of the local minima approach zero. As the 
K-Q transition is neared, = 0.57, the shape of the curve changes distinctly as the local 



minima and maxima converge. For further increases in g. Fig. 12(c) shows that the peaks 



become sharper and narrower, though the fundamental shape does not change. 

The radial distribution function (72 (r) generally shows that translational order is main- 
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FIG. 12: (Color online) Angular distribution function f{9) for various (p. 
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FIG. 13: (Color online) Radial distribution function for various (p. 

tained thorough the K and Q phases. The crystal-hke characteristics, peaks in g2{f)i are 
evident but smear out gradually as the density is reduced. Fig. [13] shows g2{r) for various q 
and several densities. While the peaks in the liquid phase, Fig. 13(c) shows that for = 0.54, 



decay quickly to unity, the peaks in the K and Q phases do not appear to decay rapidly. 
Evident in 5'2(r) for q = 2.5 and 3.5 is a split first peak that represents the subtle difference 
in pair distances of the first and second neighbors associated with the corresponding Ci lat- 
tices. In addition, the particle distribution functions, in which the particle coordinates are 
projected onto a line perpendicular to the crystal planes, are completely periodic throughout 
the K and Q phases. This suggest that lower-dimensional translational order, such as that 
in columnar phases, was not present in either the K or Q phases. 
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VI. DISCUSSION 

In this paper, we have determined the phase behavior of a general class of hard convex 
particles with shapes between a sphere and cube. We found that the degree of curvature 
plays a significant role in the cubatic ordering and phase transitions. Despite seemingly 
similar visual appearance among particle shapes, subtle differences in the curvature can 
yield large changes in the EOS. In the liquid phase, edges and corners are the dominant 
features that contribute to excluded volumes, while in the crystal phase, they can brace 
against rotations. 

Although MD has some advantages over Monte Carlo techniques, including the ability 
to simulate cooperative behavior more efficiently, there are always questions concerning 
kinetic effects and boundary effects. We have performed comprehensive tests to address 
and minimize these issues. The use of a deformable box helped to alleviate the possibility 
of anisotropic stresses and to reduce boundary effects. In each case, we verified that the 
pressure tensor remained isotropic throughout the simulations. The pressure tensor showed 
slight anisotropy immediately before the melting transition, likely due to the large stresses 
that build up and induce the transition. In all other cases, including K-Q transitions, the 
pressure tensor remained isotropic. We tested several system sizes within the limits of our 
computational capabilities and we found that the curvature of the EOS remains unchanged 
as the system size was increased. Smaller systems produced larger fluctuations and slight 
variations in the density of the melting transition compared with larger systems. 

In addition, we tested systems in what we knew to be poor initial conditions. For example, 
for q = 2.0, we initialized the system in three distinct crystals - face-centered cubic (FCC), 
simple cubic (SC), and the Ci crystal. The FCC and SC cases nearly converged to the 
pressure equation of state associated with the Ci crystal. This supports the notion that 
we have used simulation conditions that can approximate the equilibrium phase behavior. 
The internal stresses of the FCC and SC simulations were relieved when the system shifted 
toward the Ci crystal. 

As with many simulation methodologies, kinetic effects are possible due the inability to 
simulate long-time behavior. We tested several growth rates and found little variation among 
the EOS's generated for rates of I7I < 10^^. Long, constant-density trajectories revealed 
slight variation in the densities associated with freezing and melting transitions, though not 
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substantial enough to warrant concern. The curvature of the EOS in the K and Q phases 
were insensitive to the growth rates I7I < 10^^. 

The primary extension of this work includes the refinement of the boundaries in the 
phase diagram. As we have shown in Fig. [HI the DTS algorithm has shown first-order phase 
transitions quite well. However, detailed free energy calculations are needed in order to 
provide precise phase boundaries and identify higher-order phase transitions, if they ex- 
ist. In addition, we have yet to explore the dynamic behavior in hard superball systems. 
The diffusion coefficient may exhibit discontinuous jumps when crossing between K and Q 
phases. Additionally, rotational degrees of freedom play an important role in the glassy 
phase. Understanding the extent of curvature of corners and edges on the jamming char- 



acteristics has undergone initial exploration 26|. It may also be possible to manufacture 



colloidal particles with such controlled shape via photolithography or other synthetic tech- 
niques. Testing these systems for certain technologically relevant properties including wave 
scattering characteristics and rheology may reveal unusual behavior. 

While we have used a particle-growth algorithm to understand the phase behavior of hard 
particles, particle-growth algorithms are often used to search for optimal particle packings. 
In our experience, allowing particles in a dense liquid to grow slowly generally produces a 
partially crystalline system. Our results identify one of the primary challenges associated 
with searching for optimal packing arrangements since these algorithms could hardly ever 



achieve the densest state. The relaxation times are far too long for computer simu 
This highlights the need for alternative methods to find dense packings of particles 



ations. 
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28|, 



31| . Only after finding the densest packings can researchers attempt to determine the entire 



equilibrium phase behavior of hard particles by applying particle contraction. 



The application of overlap potentia 
opment of efficient MD algorithms 41 



s to generalized convex particles [43| and the devel- 



42| has made available the opportunity to explore 
more deeply how shape influences phase behavior. Along these lines, planned future work 
includes determining the onset of nematic, smectic, and possibly parquet phases of elongated 
superballs. This particular perturbation would allow one to explore the continuous evolution 
from ellipsoids to tetragonal parallelepipeds, which exhibit various types of liquid crystals 
for certain aspect ratios. With the tools available, one can determine, for example, where 
the crossover point is for the appearance of a parquet phase in a system of elongated super- 
balls. Additionally, a study of parallel hard superellipsoids, a perturbation from spheres or 
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ellipsoids to cylinders, showed the onset of a smectic phase J40|. With the addition of rota- 
tional degrees of freedom, this system that would presumably exhibit a cubatic or parquet 
phase depending on the deformation parameter. There are seemingly endless possibilities 
for hard-particle shapes and perhaps a mathematical treatment of generalized particles with 
"super exponents" may be enlightening. 
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Appendix A: Phase Behavior Effects Induced by System Geometry 

Here we report on the importance of utilizing a deformable boundary in simulations of 
anisotropic particles. Because superballs are centrally symmetric objects, we expected that 
the pressure tensor would remain isotropic throughout simulations. The EOS for the liquid 
does not depend on the geometry of simulation box. The pressure tensor for the liquid is 
isotropic since particles rotations are random. 

However, in high-density phases, we find that anisotropic stresses can build up when the 
boundaries are fixed, even in the case where q is near unity. Fig. [T3] compares the equations 
of state generated using a fixed box and a deforming box for g = 1.3 and 2.0. For small g, 
the equations of state differ significantly in their curvature. 

For q > 1.4, a fixed simulation box not only affects the curvature but also introduces a 
series of apparent first-order phase transitions. As seen in Fig. [TUfor q = 2.0, the system 
with a fixed simulation box has three apparent first order phase transitions. These dis- 
continuities along the branch are highly reproducible. Each apparent phase transition can 
be attributed to rotations of particles about certain axes. The highest- density transition 
corresponds to rotations about a single particle axis. At this density, there is enough aver- 
age spacing between particles in the lattice to allow for rotations about a single axis. The 
middle transition corresponds to rotations about two particle axes while the lowest-density 
transition on the crystal branch corresponds to free rotations about three particle axes. In 
general, the melting densities in the deforming box and fixed box differed slightly, with 
the deforming box having a lower melting density. In all cases, melting appeared to be a 
first-order phase transition. 

Observing the elements of the pressure tensor as a function of density shows that 
anisotropic stresses build up as the system approaches each apparent phase transtion. With 
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FIG. 14: (Color online) Comparison of equation of states obtained using fixed and deformable 
system boundaries for (left) q = 1.3 and (right) q = 2.0. For small q, the curvature is sharper 
with a deforming boundary than with a fixed boundary. For larger g, a series of apparent phase 
transitions are induced due to the ability of particles to rotate about certain axes. 

deforming box simulations using the Parrinello-Rahman-like algorithm allows for these in- 
ternal stress to be relaxed away quickly, as they would occur in thermodynamic systems. 
Although the Parrinello-Rahman-like algorithm does not rigorously sample an isostress en- 
semble because the lattice vectors are not directly coupled to the particle interactions, it is 
a reasonable approximation. 
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